
cd('.')
addpath('matlab functions')

D0 = readtable('NC Data/SCC_4to5.csv');

%-- Analysis--;
ra_list = [(-100:-1) (-.9:.1:.9) (1:100)];
num_ra = numel(ra_list);

clear td1 td2 fg1
for sbj = {'math','reading'}

    load(sprintf('output/main/vam_%s_matching.mat',sbj{1}),'ValueAddedEstimates','modspecs','VAmomentEstimates')
    tchcoef = modspecs.tchcoef;

    D = D0(~isnan(D0.(modspecs.outcome)),:);
    D = groupsummary(D,{'year','districtid','schoolid','teacherid'},'mean',tchcoef);
    D.Properties.VariableNames = replace(D.Properties.VariableNames,'mean_','');
    D = D(D.GroupCount>=5 & D.GroupCount<=40,:);
    [~,tchidx] = ismember(D.teacherid,[ValueAddedEstimates.teacherid]);
    vahat = bsxfun(@minus,cat(2,ValueAddedEstimates(tchidx).vahat)...
        ,VAmomentEstimates.mean);

    numTCH = size(D,1);
    VAorig = zeros([numTCH 1]);
    VAmax = struct('year',zeros([numTCH 1]),'districtid',zeros([numTCH 1]),'schoolid',zeros([numTCH 1]));
    CEVA = struct('year',zeros([numTCH num_ra]),'districtid',zeros([numTCH num_ra]),'schoolid',zeros([numTCH num_ra]));
    grpvars = {'year','districtid','schoolid'};

    for j = 1:numTCH
        VAorig(j) = D{j,tchcoef}*vahat(:,j);
        for g = 1:numel(grpvars)
            I = ismember(D(:,grpvars(1:g)),D(j,grpvars(1:g)));
            Y = D{I,tchcoef}*vahat(:,j);
            VAmax.(grpvars{g})(j) = max(Y);
            pr = ones([sum(I) 1])./sum(I);

            va_u = bsxfun(@power,exp(Y),-ra_list);
            ceva = log((pr'*va_u).^(1./-ra_list));
            ceva(ra_list==0) = pr'*Y;

            CEVA.(grpvars{g})(j,:) = ceva;
        end
    end

    [p,x] = ecdf(VAorig);
    [~,idx5] = min(abs(p-.05));
    RANKval5 = x(idx5)
    RANKorig = p(sum(bsxfun(@ge,VAorig,x'-eps),2))*100;
    RANKmax = cellfun(@(f) p(sum(bsxfun(@ge,VAmax.(f),x'-eps),2))*100,fieldnames(VAmax),'unif',0);
    RANKmax = cell2struct(RANKmax,fieldnames(VAmax));

    td1.(sbj{1}) = [mean(abs(VAorig-VAmax.schoolid)<eps);
        mean((RANKmax.schoolid-RANKorig)>5);
        mean((RANKmax.schoolid-RANKorig)>10);
        mean(RANKmax.schoolid(RANKorig<=5)>5)];

    td2.(sbj{1}) = [mean(abs(VAorig-VAmax.districtid)<eps);
        mean((RANKmax.districtid-RANKorig)>5);
        mean((RANKmax.districtid-RANKorig)>10);
        mean(RANKmax.districtid(RANKorig<=5)>5)];
   
    FIRE = cellfun(@(f) mean(CEVA.(f)<=RANKval5,1)',fieldnames(CEVA),'unif',0);
    FIRE = cell2struct(FIRE,fieldnames(CEVA));
    
    fg1.(sbj{1}) = [FIRE.year FIRE.districtid];

end

%-- RANKINGS--;
TBL = [{'Share with No Change';
        'Share with Change in Rank $>$5\%';
        'Share with Change in Rank $>$10\%';
        'Share No Longer in Bottom 5\% Rank'} ...
        num2cell([td1.math td1.reading td2.math td2.reading])];
latextbl(TBL,... 
    'tables and figures/tables/change_in_rank.tex', ...
    {'%s','%0.3f','%0.3f','%0.3f','%0.3f'});

%-- CEVA Plots--;
plotformat = {'interpreter', 'latex','FontWeight','normal','FontSize',10,'FontName','Times'};
axisformat = {'FontUnits','points','FontWeight','normal','FontSize',10,'FontName','Times'};

FOOTER = .08;
DIM = [1 2];
OUTPOS = subplot_outpos(DIM,0,FOOTER,.02);
FIGRATIO = 1.5/3;

clf
pos = get(gcf,'Position');
pos(4) = pos(3)*FIGRATIO;
set(gcf,'Position',pos)
    
i = 1;
for sbj = {'math','reading'}
    [COL,ROW] = ind2sub(flip(DIM),i);
    subplot(DIM(1),DIM(2),i)
        
    hold on
    title([upper(sbj{1}(1)) sbj{1}(2:end)],plotformat{:})
    if i ==1
    ylabel('Share with CEVA Below 5th Percentile',plotformat{:})
    end
    xlabel('$\rho$',plotformat{:})  
    set(gca,'YLim',[.00 .1],'YTick',(0.0:.01:.1),'XLim',[-100 100],'XTick',-100:25:100,'OuterPosition',OUTPOS{ROW,COL},axisformat{:})
    plot(ra_list',fg1.(sbj{1})(:,1),'--k');
    plot(ra_list',fg1.(sbj{1})(:,2),':k');
    
    i = i + 1;
end
hL = legend({'State-wide','Disrtict'},'Orientation','horizontal',plotformat{:});
set(hL,'Position',[(1-hL.Position(3))/2 (FOOTER-hL.Position(4))/2 hL.Position(3) hL.Position(4)]);

set(gcf,'PaperPositionMode', 'manual','PaperUnits','inches'...
        ,'PaperPosition',[0 0 [1 FIGRATIO]*6.5],'PaperSize',[1 FIGRATIO]*6.5)
print('-dpdf','tables and figures/graphics/va_ce.pdf','-r2000')
close all




